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ABSTRACT 

The instantaneous values of output voltages representing the wind 
velocity vector and the temperature at different elevations of the 250- 
foot meterological tower located at NASA Wallops Flight Center are 
provided with the three dimensional split-film TSI Model 1080 anemometer 
system. The output voltages are sampled at a rate of one every 5 milli- 
seconds, digitized and stored on digital magnetic tapes for a time 
period of approximately 40 minutes, with the use of a specially designed 
data acquisition system. A new calibration procedure permits the conversion 
of the digital voltages to the respective values of the temperature and 
the velocity components in a Cartesian coordinate system connected with 
the TSI probe with considerable accuracy. 

Each data-sample is divided into 58 blocks consisting of 8192 data 
points each. Stationarity of the time histories is checked by inspect- 
ing the variation of the block-means and with the use of a statistical 
test. The velocity components are transformed into the mean wind 
oriented coordinate system so that the longitudinal, lateral and verti- 
cal w'nd components are obtained. Mean values, variances and covariances 
in each data-block are first calculated and then averaged to produce 
the respective sample-means, variances and covariances. 

Power, cross, coincidence and quadrature spectra of the wind com- 
ponents and the temperature are obtained with the use of the fast 
Fourier transform. The cosine taper data window and ensemble and fre- 
quency smoothing techniques are used to provide smooth estimates of the 
spectral functions. 
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CHAPTER I 


BASIC CHARACTERISTICS OF THE ATMOSPHERIC BOUNDARY LAYER 
1.1 Description of the Lower Atmosphere 

The lower part of the atmosphere of the earth forms a boundary 
layer with a thickness of several hundred meters. The flow in the 
atmospheric boundary layer is typically turbulent, as a result of 
the combined action of shear stresses, Coriolis forces and buoyancy 
forces. The knowledge of the internal structure of the atmospheric 
boundary layer is of great importance for a large number of modern 
human activities, such as weather prediction, building and bridge 
construction, launching of aircrafts and missiles, telecommunications 
and air pollution control. 

The atmospheric boundary layer can conveniently be divided in 
smaller regions with their own particular properties. 

1. The "molecular boundary layer" is located closest to the 
ground, with a thickness of the order of one millimeter. In this 
region, molecular properties have a strong influence on the flow. 

2. The "surface boundary layer" extends up to a height between 
20 and 60 meters, depending on the conditions. It is defined as the 
region where the shear stress does not vary significantly with height, 
so that it can be considered as constant. In the surface boundary 
layer, it is mainly the surface rough. ess and the density stratifica- 
tion that determine the flow, while the Coriolis forces are negligible. 
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3. The upper part of the atmospheric boundary layer, often called 
the "Ekman layer," extends usually up to approximately 500 meters. In 
this layer, the surface conditions are of minor importance while the 
Coriolis forces cannot anymore be neglected. 

Above the atmospheric layer, the "free atmosphere" is located, 
where the air flow can be considered as inviscid and is a result of 
ine tial, Coriolis and pressure gradient forces. 


1.2 The Governing Equations 


The air in the surface boundary layer and the lower part of the 

Ekman layer can be considered as a compressible Newtonian fluid in a 

uniform gravitational nonrotating field. Let p Q , T^ and be the 

undisturbed pressure, temperature and density of the air at a certain 

height and p' t 1' and p' the respective instantaneous deviations of 

pressure, temperature and density, which can be reasonably assumed 

small, compared to p Q , T 0 and p Q . Let also U^, where i * 1,2,3, be 

the air velocity components in a Cartesian system x x x , where the 

1 2 3 

x axis is vertical. Since the flow is turbulent, it is convenient to 

3 

consider all quantities as consisting of a mean component and a fluctu- 


ating component, as follows 


u i - u i + u t 


T' * T' + 0 


p' - P' + P' 


P ■ P + P 


( 1 . 2 . 1 ) 


For short time periods, usually not exceeding one hour, the air 
flow can also be considered as steady and stationary. Steadiness of 
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Che flow implies that all mean quantities are independent of time, while 
stationarity of the flow implies that all statistical quantities are 
independent of time. 

Making use of all the previously discussed assumptions, one can 
write the general equations for the mean motion of the air in the 
surface boundary layer and the lower part of the Ekman layer in the 
following form. 

Equation of continuity 


U. . = 0 
i,i 


Momentum equation 


p U, U. . = -p , + y U. . . - u.u. . 
° j i,j , i i.jj o i 


( 1 . 2 . 2 ) 


(1.2.3) 


Energy equation 


cp r.U,=kf ..-cp 6 ,u 
p o i »ii p o >i i 


. (1.2.4) 


The dynamic viscosity, y, the thermal conductivity, k, and the 
specific heat coefficient under constant pressure, Cp, of the air are 
assumed to be constant throughout the entire region considered. 

The above equations of mean motion for turbulent flow may acquire 
a form similar to the equations of motion for laminar flow with the 
introduction of a total stress tensor, t^, and a total heat flux 
vector, H , defined respectively as 


T ij = ~ P " 5 ij + y(U i,j + U j,i ) " p o u i u j (1.2.5) 

and 


H = -k T ' , + c p 0u . (1.2.6) 

j >J P ° j 


Compared to the respective quantities for laminar flow, the above 
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introduced stress tensor and heat flux vector contain the additional 
stresses -p Q u^Uj (called Reynolds stresses) and the additional heat 
flux components c pP 0 6 U . respectively. 

The exact solution of the equations of motion in the atmospheric 
boundary layer, due to their complicated form, presents enormous diffi- 
culties. It has also been recognized that the usual mathematical analy- 
sis does not possess the capability for a clear description of the tur- 
bulent flow and, therefore, statistical and spectral analysis have to 
be used. Statistical and spectral properties have the additional advan- 
tage of being experimentally reproducible, so that an experimental veri- 
fication of related theoretical achievements is possible. 

1 . 3 The Measurement of Atmospheric Turbulence 

The properties of the lower part of the atmospheric boundary layer, 
up to a height of approximately 150 meters, are usually measured with 
instruments placed on meteorological towers or other tall structures, 
while for higher altitudes meteorological baloons and low-flying air- 
planes can be used. During the last five years, a considerable amount 
of atmospheric measurements, mostly from meteorological towers, has been 
gathered and several important conclusions have been derived from these 
experiments. A discussion of some of these tower measurements can be 
found in the references 6, 7, 17 and 26. However, these research pro- 
grams have been limited in scope in the sense that none of them provided 
a complete and multi-point set of measurements at reasonably high fre- 
quencies, Either the vertical velocity component or the fluctuating 
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temperature are missing or measurements are made with instruments whose 
frequency response is limited. 

Therefore, it is the objective of this research program to make a 
complete and extensive set of measurements of the three velocity compo- 
nents and the temperature in the lower atmosphere. This will allow the 
determination of the various statistical quantities of the velocity and 
temperature fluctuations, including the vertical momentum flux and ver- 
tical heat transfer, which are of main importance. 

In order to achieve the above requirements, measurements of the 
wind velocity components and the temperature were made at different 
elevations of the 250-foot meteorological tower located at NASA Wallops 
Flight Center. As measuring instrument, the TSI three-dimensional split- 
film anemometer was selected, which provides a satisfactory accuracy and 
a high frequency response for the velocity and temperature measurements. 

The data obtained from the measuring instruments could be processed 
either in analog or in digital form. The recent development of high-speed 
digital computers and the invention of effective digital data-processing 
techniques permitted the fast and efficient calculation of the desired 
statistical quantities from digital records of data. Therefore, a data- 
acquisition system was designed, by the use of which the analog signals 
produced by the anemometers were sampled, digitized and stored on digital 
magnetic tapes. These magnetic tapes were further processed in digital 
computers. Computer programs, based on the latest developments of the 
digital data-analysis, were used to provide all desired statistical quan- 
tities, such as mean values, variances and covariances, turbulence inten- 
sities and power and cross-spectral density functions. 
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CHAPTER II 


INSTRUMENTATION AND CALIBRATION PROCEDURES 


2.1 The TSI 1080 Three-Dimensional Split-Film Anemometer System 


The TSI (Thermo-Systems Inc.) Model 1080 Anemometer provides volt- 
age measurements, from which the instantaneous velocity vector and the 
temperature at a certain point in moving air can be calculated. The 
TSI 1080 system consists of the following components: 

1. The Probe (Model 1296E), consisting of three sensor rods for 
a total of six hot film anemometers and one thermocouple. The probe 
will be described in detail in section 2.2. 

2. One Control Circuit (Model 1053BP6T) for each probe, controlling 
the six films and the amplifier of the thermocouple. The front panel 

of the control circuit contains the pin jacks for the bridge-voltage 
outputs, the 0-5 volt outputs and the thermocouple output. In addition, 
it contains the balance potentiometers to ensure same temperature of 
adjacent films and the zero and gain adjustments for the thermocouple 
amplifier and the 0-5 volt outputs. 

3. One Voltage Regulator (Model 1051F), supplying regulated +19 
volts and -19 volts to six control circuits. 

4. One Chassis (Model 1045G6) for six control circuits. 

5. One "Rough" Power Supply (Model 1Q51FTR), with a capability 
to supply 8 regulators with ±24 volts DC and 110 volts AC. 

For further information on the anemometer and its operation, 


l 


i 


i 


1 



1 

i 


J 


1 


I 


7 


one is referred to the Operating and Service Manual for Model 1080 
Total Vector Anemometer by Thermo Systems Inc. 

2.2 The Probe 

A schematic drawing of the TSI Model 1296E Probe is provided in 
figure 1. Its main parts are the three sensors, the thermocouple, the 
supporting structure, the protective shield, the pneumatic cylinder 
and valve, the mounting flange, the electrical leads and the cable 
connector. 

The three sensor rods are mounted mutually perpendicular to form 
a Cartesian coordinate system as shown in figure 2. A thin platinum 
film of approximately 1000 angstrom thickness and 0.080 inch long is 
placed on a 0.2 inch long, 0.006 inch diameter cylindrical quartz rod. 

A thin coat of quartz over the film is used to provide environmental 
protection. The platinum film on each rod consists of two segments, 
separated from each other by two longitudinal splits 180 degrees apart. 
The orientation of the splits on each sensor rod is such that it allows 
the determination of the octant in which the instantaneous velocity 
vector is located. For minimum support interference the three sensor 
rods are mounted at only one end to a long and thin supporting structure. 

A copper-constantan thermocouple is mounted close to the sensors 
and measures ambient temperature with a maximum frequency response of 
25 hertz. 

The three sensors and the thermocouple are protected, while not 
in use, by an aluminum shield which can be placed over the sensors 
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and retracted pneumatically. Due to difficulties in supplying proper 
air pressure at the sonic orifices when the probes are mounted on a tall 
meteorological tower, the sensors are balanced at the no-flow condition. 
A weak air stream is continuously supplied through the shield, in order 
to keep the sensors free of contamination while not in operation. 

Six TSI 1296E probes were mounted on the NASA Wallops Island 
250-foot meteorological tower at the 30, 50, 100, 150, 200 and 250- 
foot levels. All probes are mounted in a vertical plane which also 
contains the axis of the tower. The maximum horizontally-measured 
distance between the two probes at the 30 and 250-foot levels is 
approximately 13 feet. Each probe is mounted with its axis in a hori- 
zontal direction, on a rotor that is controlled from the instrumentation 
trailer and can be rotated about a vertical axis to obtain any posi- 
tion between 0 and 360 degrees, with an accuracy of about i 1 degree. 

The rotor is supported by an 8 feet long boom, mounted perpendicular 
to the south face of the 250-foot meteorological tower. 

2. 3 The Data Acquisition System 

The Data Acquisition System was designed to provide simultaneous 
digitized data from one to thirteen TSI 1080 anemometer systems, sam- 
pled at a rate of 200 samples per second, for a time duration of ap- 
proximately A0 minutes. All necessary equipment and controls are 
fixed in an instrumentation trailer placed under the tower. The sys- 
tem consists of two main parts with separated operation? a) the 
mult iplex Ing and analog recording system, shown in figure A and b) the 
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demultiplexing) digitizing and digital recording system, shown in 
figure 5. 

Each anemometer provides seven analog voltages: one from the 

thermocouple and six from the hot films of each probe. All voltages 
are in the range of 0-5 volts. Each voltage is frequency modulated 
by a voltage-controlled oscillator (Data-Control Systems Inc., Model 
G0V-5C), each with a different center frequency. For the outputs 
of each anemometer there is one set of voltage-controlled oscillators 
with center frequencies of o, 12, 16, 20, 24, 28, and 32 kilohertz 
respectively. All voltage-controlled oscillators operate at a devia- 
tion from their center frequency of ±1 kilohertz for a maximum voltage 
input ±5 volts. The seven frequency modulated signals together with 
a reference signal of 100 kilohertz produced by a reference oscillator 
are summed by a Summing Amplifier (Data-Control Systems Inc., Model 
GAS-5), to produce one single multiplexed signal of approximately 3 
volts peak-to-peak. 

The six multiplexed signals from the six anemometers are simul- 
taneously recorded on six separate channels of a 14-channel analog 
tape recorder (Bell and Howell, Model VR3700B) . The time-of-day, 
provided by a time-of-day Generator (Hermes Electronics Co.) in NASA 
fast serial time code and monitored by a remote display (Itek Electro- 
Products Co., Model 2124C), is recorded on the fourteenth channel of 
the analog tape recorder. 

The demultiplexing, digitizing and digital recording system 
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tape at a time. The analog tape Is played back and the multiplexed 
signal from the desired channel Is demultiplexed to its seven components 
by passing through seven discriminators (Data-Control Systems, Inc., 

Model GFD-16). The center frequency of each discriminator corresponds 
to the center frequency of one of the voltage controlled oscillators 
and each discriminator produces an output voltage of ±5 volts. A 
reference discriminator senses any deviation in the reference fre- 
quency of 100 kilohertz caused by fluctuations in the tape speed of 
the analog tape recorder, and corrects the output of the other seven 
discriminators with a feedback procedure. 

The six outputs corresponding to the six anemometer bridge 
voltages are each passed through a low-pass filter (VTC Model LMI- 
100), which removes any frequency component over 100 hertz. This 
is necessary in order to avoid aliasing (distortion) of the velocity 
spectra at frequencies below 100 hertz. The thermocouple output 
need not be filtered since it does not contain frequencies higher than 
25 hertz. 

The seven analog voltages are fed to the multiplexing analog-to- 
digital converter (DEC Model AD01-D). It converts an analog signal 
to an 11-bit (plus sign) digital word with a time of 5 micro-seconds. 

The total scan, settling and conversion time for one sample of seven 
voltages is 50 microseconds, which is approximately 1% of the time 
interval between two samples. 

The recorded llme-of-day Js converted to a 32-bit binary coded , 

i 

decimal signal by u serial-to-paral lei converter (Hermes Electronics ) 
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Co.). This converter Is Interfaced with the computer by a Digital 
Equipment Corporation custom-designed interface. 

The Computer (DEC Model PDP 11/20) controls the multiplexing 
analog- to-digital converter and the digital tape recorder (DEC Model 
TR06-FB) through a tape transport controller (DEC Model TB68-C). The 
computer is programmable in DEC PAL-11 assembly language and has an 
8K memory of 16-bit-words. Input and output to the computer is ob- 
tained through a teletypewriter. 

The conversion starts at a time-of-day prescrioed by the opera- 
tor to assure correlated data for the different anemometers. Then 
the analog-tordigital converter performs successive scans and conver- 
sion of the seven analog voltages into two-byte words at a rate of 
one scan every 5 milliseconds. These words are stored in one of the 
four buffers of the computer. When the buffer, having a capacity 
of 209 scans or 2926 bytes, is filled with data, it is dumped on a 
9-track digital magnetic tape. A total .'f 2298 blocks or 480,282 data 
points make up one sample of data, extended over a time period of 
approximately 40 minutes. 

2.4 Conversion of the Output Signals into the Actual Quantities 

A detailed discussion of the theory of operation of hot-wire and 
hot-film anemometers can be found in the references 15 and 22. This 
report is limited to a discussion of the empirical expressions obtained 
through the calibration of the probes from which the three velocity 
components in the sensor-oriented coordinate system can be calculated. 
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Each one of the six films of the probe is connected in a bridge 
circuit as shown In figure 6 and operates at a constant temperature. 
Special care has been taken so that both films of each sensor operate 
at the same temperature, In order to prevent heat transfer from one 
film to the other. 

The bridge voltage is amplified in a feedback loop to ascertain 

thermal inertia compensation and is measured to provide the value E . 

ai 

where the subscript a • A,B,C stands for the ser e. rod and the sub- 
script i ■ 1,2 stands for the appropriate film 0 ilgure 2 for the 
exact location and orientation of all six films). The 0-5 volt out- 
put circuit is set up so that the bridge voltage for each film can 

be obtained from the corresponding 0-5 output voltage, V .with 

ai 

the following expression 


E ai " 2 * V a i + 3 * a * A » B ’ C * i - 1.2 


(2.4.1) 


The h«at convected from the electric heated element to the ambient 
air is proportional to the square of the voltage drop in this element 
and can be approximated by the following expression 

Q a- C V K a.- E a5 + C W E „i' « ' A ' 8 ' C • 


where K and K are calibration constants for each film and CF 
al a2 a) 

and CF are correction factors to account for cable and ambient tempera- 
ture changes. According to reference 22, the expression for these 
correction factors is 
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The zero-superscribed quantities refer to the TSI calibration 

conditions. The change of the cable resistance due to a change in 

cable temperature from the TSI ambient temperature, T °, is equal <n 

& 


AR 


-R 0 a °(T - T °) 


(2.4.4) 


"ai L ai ^ai w ai 
where the temperature coefficient of resistivity for copper, refer- 
enced to the temperature T a °, can be calculated as 


1 


cu 390.1 + T 


(2.4.5) 


m 


It R is the resistance of the cable measured at a tempera- 
ai 

ture, T c m , then the cablt resistance at TSI ambient temperature can 
be calculated from the following expression 


R ° - R m 
C ai c ai 


1 + 


390.1 + T m 8 
c 


(T 0 - T B ) 
c 


(2.4.6) 


The hot-film resistance, R , at TSI conditions can be calculated 

f ai CQld 

from the measured cold-film resistance, R^ , as 

f ai 


i 0 - Rjoldjj + a cold (T o _ T cold 

ai f ai ai f a 


)] 


(2.4.7) 


where a£ 0 ^ is the measured temperature coefficient of resistivity o£ 
ai 

the film. 

The hot-film resistance at operation conditions in the wind 
tunnel or in the atmosphere can be obtained from the hot-film re- 
sistance at TSI conditions according to the expression 


R 


+ Ar 


ai ^ai '"ai 
This is due to the fact that the total resistance in the 


(2.4.6) 
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bridge arm containing the cable and the film must be constant. There- 
fore, a change in cable resistance causes a change in hot-film resis- 
tance and, as a result, changes the operation temperature of the film. 
The new operation temperature of the film can be calculated from the 
following expression 


A R 


= T 0 + 


-ai 


f ai f a a cold R« ld 


(2.4.9) 


f ai f ai 


In order to be able to use the same empirical expressions when 
the density of the air changes, all velocities are calculated for 
standard pressure and temperature of 

p, ■ L4.7 psia * 29.92 in Hg 
T g - 530° R 

Since the hot-film sensors measure the mass flux, pU, where p 
is the air density, the standard velocity U , can be obtained with 
use of the equation of state as follows 


T 

U = U * . (2.4.10) 

S Ps T a 

When the hot film is heated in flowing air, the heat convected 
from the film to the air depends on the velocity magnitude as well 
as the angle of attack, <J>. The effective cooling velocity, U £ , is 
defined as the normal velocity that would produce the same amount 
of cooling from the film, as produced by the actual velocity, U, 
at a particular angle of attack, <J>. The standard effective-cool- 
ing velocity, U eS , can be obtained from the heat transfer,. Q n , 
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through an expression of the following form 


a D .u n 

T ° - T 0 
f o 3 


(2. A. 11) 


where the constants D and n are to be obtained through calibration 
of the probes. T.n order to obtain identical calibration curves (heat 
transfer versus standard effective-cooling velocities) for all three 
sensors on the same probe, the calibration constants K ^ need to be 
adjusted. 

According to reference 8, the standard effective-cooling velocities 
are related to the magnitude of the standard velocity, U g> and the 
angle of attack, <+>, by the expression 


U 2 * U s 2 (cos 2 <ji + k sin 2 4> a ) , a * A,B,C , (2.4.12) 

a a 

where the constants k^, kg, k^ are of the order of 0.1 and have to be 
determined through calibration. The three sensor-yaw angles <^, $>g 
and have the following trigonometric relationships 

sin 2 4> A + sin 2 <J>g + sin 2 ^ = 1 

and (2.4.13) 

cos 2 4> A + cos 2 $g + cos 2 4>£ * 2 

Assuming that identical values of k^ exist for all three 
sensors and using the above trigonometric relationships, one obtains 
the magnitude of the standard velocity in terms of the standard 
effective-cooling velocities of the three sensors, as follows 


5 



mm-m 


i 


j 


1 


i 



I 


u 2 + u 2 + u 2 

eS A eS B es C 
k„„ 


(2.4.14) 


where the constant k av is of the order of 2.1. Now the sensor-yaw 
angles can be calculated from 

r r u 


|<f> a l “ arcsin 


1 - C 


1 - k 


(2.4.15) 


where the coefficients C A> C g , C c are the reciprocal of the respec- 
tive maximum values of U /U at zero sensor-yaw angle. Ideally, 

6S S 

a 

the values of these coefficients should be equal to one, but due to 
experimental errors and errors in curve fitting they may deviate a 
small amount from unity. 

The sign rf each angle can be determined by comparing the ratio 
of the voltages of the two films of some other sensor to the ratio, 

R, of the voltages of the same two films, when the velocity vector lies 
in the plane of their split. Hence, 



*A > 0 

when 

V E C2 

>R c 

9 

1 

♦b > 0 

when 

V E A2 



1 and 

! 

! 

*C > 0 

when 

E BL /E B2 

>r b 

# 


The actual velocity components in the sensor-oriented coordinate 
system will be calculated by the eauations 

U a - U sin<(> a , a =* A,B,C , (2.4.16) 

where the actual velocity magnitude, U, is calculated from the standard- 
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velocity magnitude, U g , by 


P T 

U - U * — * — 
s p T s 


(2.' .i/) 


The ambient temperature can be calculated from the output voltage 
of the thermocouple, according to the following linear expression 


C T 1* V T + C T2 


(2.4.18) 


where the values of the constants C and C are determined through 

T1 T2 

calibration of the thermocouple. It is assumed that the thermocouples 
are not affected by the heat convected from the hot films. 


2.5 Calibration of the Anemometers 

It was concluded in reference 22 that the data analysis suggested 
by TSI does not possess the required accuracy. Therefore, a new 
method based on new calibration procedures was developed. It was also 
concluded that accurate measurements with the TSI 1296E probe are only 
obtained when the mean velocities are in a direction nearly-parallel 
to the axis of the probe. This can be made possible for the tower 
measurements by rotating the probe, until it is approximately in the 
mean wind direction. If during the measurement period of 40 minutes 
the mean wind would change direction appreciably, the data would be 
discarded as being non-stationary . 

The calibration of the split-film sensors was carried out in 
the low-speed Aerolab wind tunnel, located in the Quality Veri- 
fication and Calibration Facility at NASA, Wallops Flight Center.. 
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The test section dimensions are 32 inches (height) by 45 inches 
(width) by 48 inches (length). The air flow is produced by a three- 
blade variable pitch propeller, moved by a 1160 RPM / 50 HP con- 
stant speed electric motor. Temperature and air speed in the test 
section are measured with an electronic thermometer (Atkins H-51) 
and twc electronic manometers (CGS Models 1015 and 1Q14A) respec- 
tively. Their outputs, together with the six voltages from the 
TSI anemometers are scanned, digitized and averaged for a period 
of approximately 2 seconds. The data acquisition system is con- 
trolled by a HP 9810A programmable calculator and a HP 2570A 
coupler /controller. The average values of the temperature, the 
standard velocity and the six bridge voltages from the TSI anemo- 
meter are printed on paper tape. 

The probe was mounted on a vertical support and could be 
rotated in a horizontal plane at any probe-yaw angle, 8, while 
the probe-pitch angle, a, was held equal to zero. Rotation of 
90 degrees of the probe on the support offered the possibility 
of varying the pitch angle, a, while keeping the probe-yaw angle, 

8, equal to zero. Three sets of data, each for 13 different flow 
velocities between 3.5 and 50 feet per second were first taken, 
with the probe in such a position that one sensor at a time was 
perpendicul.ar to the flow. From the geometry of the sensor array 
it can be seen that A is perpendicular to the flow when a - -35.26 
degrees and 8=0 degrees, B is perpendicular to the flow when 
a « 0 degrees and 8 = 39.25 degrees and C is perpendicular to the 
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flow when a » 0 degrees and $ ■ -39.25 degrees. Nine more sets of 
data, each for five air velocities between 6 and 35 feet per second, 
were also taken for a - 0 and B ■ 0, ±10, ±20, ±30 and ±50 degrees. 

The data required for the calculation of the temperature cor- 
rection factors for each film as given by equation 2.4.3 are ob- 
tained from the following measurements. The resistance of the cold 
films, R co ^, at ambient temperature T C °^, as well as the cable 


ai 


resistances, R , at temperature T , are measured with a preci- 
c ai c 

sion digital volt-ohm-meter (Fluke Model 8200A) . The temperature 

coefficient of resistivity of the films, a co ^, is calculated 

^ai 

from measuremei. s in the thermal-chamber, also located in the 
Quality Verification and Calibration Facility. 

The different constants that are necessary for the calcula- 
tion of the velocity components will now have to be determined from 
the above measurements. First, the rate at which heat is convected 
from each sensor when it is perpendicular to the flow is calculated, 
by making use of the equation 2.4.2 and using values suggested 
by TSI as a first approximation for the constants K^. Then the 
heat fluxes for the three sensors are plotted versus the standard 
velocities and the values of K ^ are adjusted, so that the data 
collapse to one single curve. 

The next step is to determine the proper coefficient D and 
exponent n in the expression 2.4.11, which will fit the experi- 
mental data. It has been found that not one single set of values 
for D and n will describe the calibration curve adequately. The 
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calibration curve has to be divided Into three ranges with different 
values of D and n in each range. 

Once the heat flux versus effective cooling-velocity relation 
has been established, the standard cooling velccities, U gg , can be 
calculated for different sensor-yaw angles. Accepting that an 
expression of the form of 2.4.14 exists, an average value of the 
coefficient k flV can be obtained, for which the probe is expected to 
operate most of the time. The value of the standard velocity magni- 
tude, V g , is obtained from the measurement of the dynamic pressure 
by a Pitot-static tube in the wind tunnel. 

The final step is to determine the constants for the expression 
2.4.15, which yield the magnitudes of the sensor-yaw angles. These 
constants for sensors B and C are obtained by fitting a curve of the 
form 2.4.15 to the calibration data, obtained by plotting of the 
velocity ratios U /U versus the respective sensor-yaw angles. In 
order to obtain the yaw angle for sensor A, use is made of the trigono- 
metric relationship 2.4.13. 

The copper-constantan thermocouple is independently calibrated. 
First the zero and gain potentiometers were set up approximately, so 
that an output-voltage range of 0 to 5 volts DC corresponds to a 
temperature range of 0 to 200 degrees Fahrenheit. Then the thermo- 
couple was placed in ther thermal-chamber and measurements of the out- 
put voltage were taken for temperatures in the range from 20 to 100 
degrees Fahrenheit. Fitting the so obtained data to a straight line 
permits the calculation of the values of the constants in the expression 


2.4.18. 
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As an example of the above calibration procedure and its accuracy, 
the calibration constants and the observed calibration errors are pro- 
vided for the TSI #1193 Anemometer, which was placed in the 100-foot 
level of the meteorological tower. The experimental values of U eg /U s 
and the semi-empirical relation of the form of equation 2.4.15 are 
plotted versus the sensor-yaw angle, #, in figure 7 for sensor B and in 
figure 8 for sensor C. The calibration curve for the heat flux versus 
the standard effective cooling velocity for all three sensors is shown 
in figure 9. All necessary calibration constants are provided in 
table I. Finally, the observed calibration error in the calculation 
of the velocity components versus sensor-yaw angle is shown in figure 
10 and the error in the calculation of the velocity magnitude versus the 
probe-yaw angle is shown in figure 11. 
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CHAPTER III 

STATISTICAL ANALYSIS OF THE DATA 
3. 1 Limitations of Discrete Time Histories 

Consider a physical quantity x(t), which shows a random variation 
with time. Every record, x (t), of the output of an instrument measur- 
ing the quantity x(t) over a finite time interval is called a time 
history. Every record of x(t) extending over an infinite time interval 
is called a sample function. The ensemble of all possible sample func- 
tions, that can be recorded if the measurement of x(t) is repeated an 
infinite number of times, is called time series (or random process 
or stochastic process) and is denoted as {x k (t)}. The assumption of 
ergodicity permits the derivation of desired information about the 
entire time series from the analysis of a single arbitrary time history 

In many cases, the continuous output of the measuring instruments 
is digitized and sampled at some appropriate uniform time increment, At 
over a total time interval, T. Then the time series which is being 
measurfd is represented by the discrete time history, where by 

definition 

x = x(ivit) for i «* 1,2,. ..,N (3.1.1) 

i 

and the total record time, T, is equal to the time increment, At, times 
the number of data points in the sample, N. 

The analysis of discrete time histories is subjected to the fol- 
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lowing constraints: 

1) Sampling of the time aeries. 

11) Finite record length, 
ill) Discretization of data. 

iv) Restricted range of data, due to the limited response of the 
transducer used. 

It is required that the effects of these constraints are confined 
by the use of appropriate techniques. 


3.2 Basic Statistical Definitions 


Giver any sample function of time, x(t), its mean value, p , is 
defined as 


u x * litn Tr f X(t)dt 
X T-** ZT „T 


(3.2.1) 


while its mean square value, 4^, is defined as 


4^ ■ lim i_ f x 2 (t)dt 
T-*° 2T _ T 


(3.2.2) 


The average of the squared differences of each value of x(t) from 
the mean value is called the variance, c 2 , defined as 


3 2 ■ lim / t x (t) - V J 2 dt 

T-+00 2T -T * 


(3.2.3) 


The square root of the variance, o x , is called standard deviation. 
Given two sample functions, x(t) and y(t), their covariance, o 2 , 

xy 


is defined as 


°xv ■ lim if ; tx(t) ’ ' ^J dt 

7 T-kc -T 7 


(3.2.4) 
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The dependence of future values of the considered sample functions 
x(t) or y(t) upon the present values of x(t) is described by the auto- 
correlation function, R x (r), and the cross-correlation function, R xy (i), 
which are respectively defined as 

T 


and 


Rx(t) ■ lim -- / x(t) x(t + T)dt 

T-*» -T 


R^Xt) - lim i— / x (t) y(t + T)dt 
3 2T _t> 


(3.2.5) 


(3.2...) 


When a discrete time history, x^ where i ■ 1,2,...,N is given, 
its mean value, x, is defined as 

N 


- 1 - 
11 *5 E x i 

i-1 


(3.2.7) 


Then any instantaneous value, x^, can be considered as the com- 
bination of a static or mean component, x, and a dynamic or fluctuating 
component, x^» so that 

x t - x + x i . (3.2.8) 

As a result of the definition of the mean value, it comes out that 
the mean value of the fluctuating component is equal to 2 ero, i.e. 

X - 0. 

In many cases it is necessary to remove the mean component from 
the discrete time history and therefore generate a new mean-free time 
history, which contains only the fluctuating component. This is ob- 
tained by subtracting the mean value, x, from each instantaneous value, 
x^, of the original time history. 

If x^, for i ■ 1,2,. ..,N, is such a mean-free time history, its 
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variance, x 7 , Is given by the expression 

_ N 

it 7 ■ ^ [ x 2 . (3.2.9) 

N i-1 1 

If x and v are both mean-free time histories, their covariance, 
i 'i 

xy, is defined as 

- 1 N 

xy - - l x.y . (3.2.10) 

N i-1 1 

The existence of slowly varying mean values in the time histories 
results in great distortions of their spectral estimates. Therefore, 
it is good practice to remove the means from the data before proceeding 
with the spectral analysis. The removal of the means proved necessary 
for time histories related to atmospheric measurements, such as air 
velocity, temperature and pressure. 

3.3 Partition of Each Sample of Data into Blocks 

As described in Chapter II, each sample of data from each TSI 1080 
probe consists of four discrete time histories, namely the three velo- 
city components in the sensor oriented coordinate system and the temper- 
ature, sampled at a rate of 200 points per second, for a time period 
of approximately 40 minutes. The total number of data points in each 
time history is approximately 48C 3 000. 

The simultaneous spectral analysis of this number of data points 
is far beyond the capacity of any available computer. It is therefore 
necessary that each time history is divided into a number of data blocks 
and that the statistical and spectral analysis is applied separately to 
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each block. As will be explained in Chapter 4.3, the most efficient 
use of the fast Fourier transform requires that the number of data 
points per block is equal to an inte^... power of 2. A: c tur considera- 
tion of the capacity of the IBM 370 computer that was used, the number 
of data points per block was chosen to be N » 2 1 3 ■ 8192. Consequently, 
the number of data blocks in each time history had to be taken as M * 58, 
so that the total number of data in each time history was reduced to 
M*N » 475,136 , corresponding to a time period of approximately 39 min- 
utes and 35 seconds. 

Assuming that the data sample is stationary, mean quantities, vari- 
ances and covariances in each block w*re first calculated and later 
all block means, variances and covariances were averaged to provide the 
sample mean values, the sample variances and the sample covariances of 
the fluctuations of the three velocity components and the temperature. 

3.4 Stationarity Tests 

As already mentioned in former chapters, the statistical and spec- 
tral analysis is based on the assumption that the time histories are 
stationary. Therefore, it is required that this assumption is carefully 
checked. This will be done in two ways: 

1. By simply inspecting the variation of the block means, the 
block standard deviations and the block probe-yaw angles. Large varia- 
tions of the mean wind, while large variations of the standard devia- 
tions indicate that the frequency structure of the turbulence is vary- 
ing with time. The magnitude of the probe-yaw angle, i.e. the angle « ' 
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between the longitudinal axis of the probe and the direction of the 
horizontal mean wind, will also determine the level of accuracy of the 
probe operation, since the most accurate data are obtained when the 
probe is directed straight into the mean w'nd direction. 

2. By detecting monoconic trends using the following nonstation- 
aritv trend test. 

If A*, i * 1,2,...,M represent any uf the statistical quantities 
calculated in the i_ data block, the total number of reverse arrange- 
ments, R, is defined as 

M-l M 

R = 1 I R,, 

i=l j-i+1 

where 

1 , if A 1 > 

0 , if A 1 £ A j 

Assuming a normal distribution for the number of reverse arrange- 
ments, there exists a certain interval [R ,R ] for the values of R, in 

1 2 

which one can be assured at a certain level of confidence, u, that 
th° data sample does not show any monotonic tiend. ^or values of R 
greater uhan R , a downward trend of the statistical quantitv tested 
is to be expected, while for values of R smaller than R , the data 
samnle shows an upward trend of the respective statistical quantity. 

The limits of the stationarity interval have been calculated in 
reference 23 for a total number of blocks, M, between 10 and 100. 

For M - 58, it was found that at a confidence level, a, equal to 90%, 

the lower limit, R , is equal to 703 and the upper limit, R , is equal 

1 2 



(3.4.1) 


(3.4.2) 
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to 949, while at a confidence level 98%, R is equal to 652 and R is 

1 2 

equal to 1000. 

The stationarity trend test Is performed for the block means of 
the. three velocity components, the magnitude of the velocity vector 
and the probe-yaw angle, as well as the block standard deviations of 
the three velocity components. 

Based on the results f the inspection and the nonstationaritv 
trend test of the block statistical quantities, a decision will he 
made at this point whether or not to continue the statistical analysis 
for each one sample of data. 

3 . 5 Coordinate Transformation into the Mean-Wind Oriented Coordinate 

System 

The TSI 1030 anemometers provide the three components U^, U R and 
U c of the wind velocity vector in the Cartesian system that is formed 
by the three sensors, as shown in figure 2. 

Nevertheless, the velocity components have to be expressed in some 
coordinate system, which is physically connected with the flow field. 

For this purpose, the mean-wind oriented coordinate system, with co- 
ordinate directions x,y and z was introduced. This is a Cartesian 
system, where the z-axis is vertically upward, the x- and y-axes are 
horizontal and the sample mean velocity vector has to lie in the x-z 
plane. Then, bv definition, the y-componcnt of the sample mean velocity 
is zero. On the other hand, since the flow in the atmospheric boundary 
layer is almost horizontal, it is p.xpected that the z-component of the 
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sample mean velocity vector will be very small compared to the x-component. 

Using geometrical relations between the two coordinate systems, 
one finds that the transformation will have the following matrix form 

U. 


(“1 


■V 

1= 

w 



Ur 


0. 57735cosB O.57735cos0-O.7O711sinB 0.57735cosp+0. 70711sinp 

0. 57735sinB 0. 57735sinB+0.70711cosB 0.57735sinB-0.70711cosB 

-0.81650 0.40824 0.40824 

(3.5.1) 

where u, v and w are the velocity components in the xyz-system and P is 
the sample probe-yaw angle, defined as the angle between the axis of the 
probe and the x-axis, as shown in figure 3. 

According to the geometry of the probe, the sample probe-yaw angle 
can be calculated from the following expression 

(3.5.2) 


B * arctan 


1.22475 


U 0 - U 
C B 


U A + U B + U C 


where U^, U^ and U are the sample mean values of the velocity components 
U^, and U^ respectively. 

3.6 Calculation of the Mean Values, Variances and Covariances 


After the transformation of the velocity components into the mean- 
wind oriented coordinate system, the time histories to be analyzed are 
the temperature, T (in degrees Fahrenheit), the longitudinal velocity 
component, u, the lateral velocity component, v, and the vertical velo- 
city component, w (all velocity components are expressed in feet per 
second). Again as before, each time history is divided into 58 blocks 
of 8192 data points each. 
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The sample-mean values are calculated in two steps. First, the 
block-mean values, 5°, are calculated from the expression 


N 


- n 1 

X “17 lX 

i=l 


N . , i 


(3.6.1) 


and then all block-means are averaged to produce the respective sample- 
mean values, x, as follows 

M 


1 - n 
x ■ « 


(3.6.2) 


where x represents any one of the quantities T, u, v or w. 

For the calculation of the variances and covariances, it is neces- 
sary to use time histories with zero mean values. Therefore, the 
block-means have to be subtracted from all values in each block and 

■■ - “rt l ■ ■ "ft 

the block-variances, x , and sample-covariances, xy , of the time 
histories represented by x and y, are calculated as 


N 


^ i z - x n ) 


and 


N 


xy 


N i^ 1 


(3.6.3) 


» Z (x . - x n )(y. - y") 
i»l 


(3.6.4) 


The sample-variances, x 7 , and sample-covariances, xy, are calculated 
as the average of the respective block-quantities, as follows 


M _ 

- I Z 
x M t x 

n=l 


(3.6.5) 


ami 


xy = 


M 


)’ xy 

n-1 


(3.6.6) 
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CHAPTER IV 

SPECTRAL CALCULATIONS 

4 . 1 Definitions 

The nature »nd properties of a physical system, whose character- 
istics are measured as functions of time, are well illustrated with 
the spectral density functions in the frequency domain. 

A function that describes the general frequency composition of the 
measured quantity, x(t), and the relative significance of every frequency 
is the power spectral density function, G x (f). There exist three ways 
by which the power spectral density function can be defined. 

1. From the direct Fourier transform of the function x(t), accord- 
ing to the expression 


G (f) = 2’lim 
x T-*» 


/ , s -12TTft 
/ x(t) *e dt 

-T 


(4.1.1) 


2. From the Fourier transform of the autocorrelation function, as 
foil ows 


G x (f) =2 / R x (t) *e i “ 7Ift dT 

—00 


(4.1.2) 


3. From the mean square value of the output of a sharp band-pass 
filter. The filter output, x(t,f,Af), consists of the portion of x(t) 
that lies in a frequency range with center frequency f and bandwidth Af , 
The power spectral density function is then defined as follows 
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G < f > “ lim 77 lim ~ f x 2 (t, 
Af-K) a L T-«° -T 


(4.1.3) 


All three definitions of the power spectral density function are 
asymptotically equivalent, when applied to stationary time series. The 
power spectral density function is always a real-valued, non-negative 
function of frequency. Its complete plot versus frequency is called 
power spectrum. Since negative frequencies do not possess any physical 
meaning, the power spectrum extends only in the positive frequency domain. 

The joint properties of two functions of time, x(t) and y(t), in 
the frequency domain, are described by the cross-spectral density func- 
tion G (f). Similar to the power spectral density function, the 
cross-spectral density function can be defined in three equivalent ways. 

1. From the Fourier transform, Y(f,T), of y(t) and the complex 
conjugate, X*(f,T), of the Fourier transform of x(t), as follow - 


G xy (f) = 2 lim X*(f,T)-Y(f,T) 
' T-*” 


, (4.1.4) 


where the finite Fourier transform, X(f,T), of the function x(t) in the 
iacerval (-T,T) is defined as 


X(f,T) = / x(t)*e ^’^ 1T ^ t dt 

-T 


(4.1.5) 


2. From the crcss-correlation function, according to the expres- 


G xy (f) - 2 f R xy (i) e" i27Tft dT 
— 00 


(4.1.6) 


3. From the product of the output of two band-pass filters, ap- 


i 


i 


1 



1 
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plied to the £unc..j.oa> x(t) and y(t), according to the expression 


G xv (f) * lim XT ' lim W f x(t,f,Af).y(t,f,Af)dt - 
y Af-0 T-*» ZT -T 


-i / x(t,f,Af)-y°(t,f,Af)dt 
-T 


(4.1.7) 


where the function y°(t,f,Af) is equal to the function y(t,f,Af) shifted 
90 degrees in phase angle. 

The cross-spectral density function is in general a complex func- 
tion of frequency, consisting of a real part, C X y(f) , called coincident 
or co-spectral density function and an imaginary part, -Qxy(f)» called 
quadrature or quad-spectral density function, i.e. 


G xy^ “ G xy^ ” iQxy(f) 


(4.1.8) 


Alternatively, the cross-spectral density function can be expressed 


in polar form as 


G xy< f > “ l G xy( f >l‘ e 


-i6™(f) 


(4.1.9) 


where the magnitude and the phase angle are related to the co- and 


quad-spectral density functions by the expressions 

IV £)I ' [c ^ (£) +Q ^y (£) ^ 


(4.1.10) 


Q (f) 

6(f) = arctan ~^L — 

^ lc xy (f) 


(4.1.11) 


A real-valued quantity, which is related to the cross-spectral 
density I'unc l ion, is the coherence function, defined as 


Y 2 U) 

I VV ' * 


l G xv (£) l 2 


G x (f)-G y (f) 


(4.1.12) 
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The coherence function is a measure of the correlation between the 
values of two functions of time and takes values between 0 and 1. When 
the coherence function is equal to zero for a certain frequency, f, the 
two functions, x(t) and y(t), are said to be incoherent at the frequency 
f. When the coherence function is equal to unity for all frequencies, 
the two functions are said to be fully coherent. 

The graphical representation of the cross-spectral density function 
is called cross- spectrum. It consists of the plot of both the magni- 
tude and the phase angle of the cross-spectral density function versus 
frequency. 

Based on the definitions mentioned above, three different methods 
for the calculation of the power- and cross-spectral density functions 
are available. The recent development of FFT (fast Fourier transform) 
techniques made the direct Fourier transform exceedingly advantageous 
for the spectral analysis of long time histories. Therefore, the 
direct Fourier transform method based on a new FFT technique, as des- 
cribed in the following sections, will be used. 

4.2 Finite Interval Effect — Cosine Taper Data Window 

Any time history is by definition extended over a finite time inter- 
val. Definitions and equations, which are proper for infinite time 
series, may be erroneous, when applied to time histories. Therefore, 
careful consideration of the finite interval effect has to be made. 

The genera] Fourier transform, X(f), of the time function x(t) is 
defined as 

\ 
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X(f) « f x(t)*e i ^ 1T ^ t dt 


(4.2.1) 


while the finite Fourier transform, X(f,T), of the time history x(t), 
where t lies in the interval [-T,T] , is defined as 


X(f,T) - / x(t)-e" 127rft dt 

-T 


(4.2.2) 


Any time history can be considered as produced by a sample function, 
x(t), multiplied by a window function, w(t), which eliminates all values 
of x(t) outside the interval [ — T , T] . The window function that permits 
all values in the interval [-T,T] to pass undeformed is called the 
rectangular or boxcar window function and is defined as 

1 | t| <_ T 

0 T <_ | t| . (4.2.3) 

With the use of the rectangular window function, the finite Fourier 
transform can be alternatively expressed as 


w(t) 


X(f,T) ■= f x(t)-w(t)-e' i27rft dt 


(4.2.4) 


According to the convolution theorem, the finite Fourier transform can 
also be considered as the convolution of the general Fourier transform, 
X (f ) , with the Fourier transform, W(f,T), of the window function, i.e. 


X(f,T) - f X(f - 0*W(C ,T)dS 


(4.2.5) 


The Fourier transform of the rectangular window function is cal- 
culated lo he 

, r v Z — 12nft. r -i2 n f t . sin 2 TTf t 

W(f ,T) ■ /w(t)* e dt » J e dt * . (4.2.6) 

-T "f 
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The rectangular window function and its Fourier transform for 
positive frequencies are shown in figure 12. 

The effect of the rectangular window can be illustrated with the 
following example. Consider the sinusoidal function x(t) » cos 2irf Q t, 
where f is a constant frequency. Its infinite Fourier transform is 

\ calculated to be 

f 

I x(f) - ms(f + f ) + 6(f - f )] , (4.2.7) 

; o 

| while its finite Fourier transform in the interval [-T,T] is calculated, 

according to equation 4.2.5, as 

’ sin 2 tt (f + f )T sin 2n(f - f )T * 

X(f,T) = h — + 2— . (4.2.8) 

n(f + f J ir(f " f ) 

L 0 0 J 

The functions X(f) and X(f,T) are plotted versus frequency in 

figure 13 for the range of positive frequencies and for the case when 
the product f^T is equal to 5. As can be seen in this figure, the 
function X(f) has a non-zero value only at the point f ■ f Q , while the 
function X(f,T) extends over the entire frequency domain and presents 
positive and negative peaks, called sidelobes. 

When x(t) is a wide-band function (i.e. having a wide frequency 
spectrum), its finite Fourier transform will consist of the superim- 
posed effects of all frequencies, according to equation 4.2.5. Al- 
though the negative sidelobes effect will now be much smaller than in 

t 

| , the narrow-band function case, a distortion in the spectral estimates 

| is to bo expected. Under certain circumstances, it is even possible that 

? 

r 

I negative values of the power spectral density function appear for 



certain frequencies, which is opposed to the physical significance of 
this function. This danger is reduced by replacing the rectangular 
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window function with some other window function, which presents lower 
sidelobes. 

From the existing number of data window functions, which have been 
suggested and used by different authors, the cosine taper data window 
function was chosen as most satisfactory for meteorological time his- 
tories. The cosine taper data window function deforms one tenth of the 
data at each end by multiplying these data by a half cosine bell and 
leaves the rest eight tenths of the data undeformed. This window func- 
tion is defined in the interval [-T,T] as follows, 

fo. . T < |t| 

0.8T <_ 1 1| <_ T 

V ) 

1 0 <. 1 1 1 < 0.8T . (4.2.9) 

. 

The Fourier transform of the cosine taper data window function is 

(25 - 12 f 2 T 2 ) [sin(2irfT) + sin(2irf0 .8T) ] 

W(f,T) . (4.2.10) 

2irf [25 - 4f 2 T 2 ] 

Both the cosine taper data window function and its Fourier trans- 
form are shown in figure 14. Compared to the rectangular window function, 
it shows lower sidelobes and, therefore, produces less leakage in the 
spectral estimates. Nevertheless, the spectral estimates have to be 
multiplied by a scale factor of 1/0.875, due to the deformation of part 
of the data. 

4.3 The Fast Fourier Transform 
4.3.1 Introduction 

Consider a real- or complex-valued discrete time sequence, x^, 
for k ■ 0,1,2, ... ,N-1 which consists of N data points sampled at a 


w(t) 


1 + cos 


10*t 

IT 


I 


1 



1 


1 


I 
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uniform time increment. At, so that the total time Interval, T, is 
equal to N*At. The finite Fourier transform of this time sequence, 
as defined by equation 4.2.2, has to be approximated by the discrete 
Fourier transform, X r , defined as follows 


N-l -i2irl 


X - ^ l x e 
r N k-0 k 


, r - 0,1, 


(4.3.1) 


This expression can be obtained from equation 4.2.2 by replacing 
kT T 

the time, t, by — , the time increment, At, by — and the frequency, f , 
by The discrete Fourier transform consists of N discrete values, 
corresponding to N consecutive frequencies, separated from each other 
with a frequency increment, Af, equal to 

The computation of the discrete Fourier transform of a complex data 
sequence of N points, according to equation 4.3.1, requires N complex 
multiplications and add ' tions . For an adequately large number of data 
points, this nuaber of operations would consume an enormous amount of 
computer time and makes equation 4.3.1 practically useless. A more 
efficient method for the computation of the discrete Fourier transform 
of long data sequences is the fast Fourier transform, requiring only 

2Nlog N complex operations, instead of N 2 . For 8192 • 2 13 data points 

2 

the FFT requires approximately 315 times less computer time than required 
by using the equation 4.3.1. 

There exist several FFT techniques, varying in their details. A 
technique, having the additional advantage of not requiring the bit- 
reversal of the final values, will be presented in the following section. 
4.3.2 Derivation of the Equations 

The highest efficiency of the FFT is obtained, when the number, N, 


i 


J 


i 


i 


! 


1 


1 


i 

i 




39 


of data points to be transformed is equal to an integer power of 2. If 
this is not the case, the number of data can either be truncated to the 
nearest lower power of 2 by omitting some of the data or reach the near- 
est higher power of 2 by adding some zero-valued data at the end of the 
sequence. 

p 

Consider a time sequence of N points, where N is equal to 2 , 
p » 2,3,..., extending over the time interval [0,T]. Letting W ■ e~* 2 
for convenience, the equation 4.3.1 can be rewritten as 


N-l 

X - i Ex.«“ , r - 0,1 N-l 

r N k-0 K 


(4.3.2) 

The integers k and r can be expressed in the binary number system as 


k-h ,-2P _1 +h ,*2 p_2 +-‘-+h *2^ *2 l +h n *2° 

P-1 P-2 2 1° 

, h A - 0,1 (4.3.3) 

and 


r - j p-1 •2 P_1 +j p _ 2 *2P”2+. . .+j ^«2 2 +j ^ •2 1 +j 0 *2° 

, - 0,1, (4.3.4) 

or, according to a different notation, as 


k - (hp-i » * • • ,h Q ) 

(4.3.5) 

and 


r * <Vl j i ,J o > 

(4.3.6) 


Following this notation, equation 4.3.2 can be written as 

T N-l 4 

X(J j ,j ) — — lx (h . , < • • *h »h Q ) W , (4.3.7) 

p 1 0 N koQ V 1 


or, furthermore, 


x ' Vi VV " m 1 


1 1 

y 


rk 

.U 


l x(h j,...,h ,h )W .(4.3.8) 
h -0 h -0 h i-0 J 

O j p A 


The fLrst step of the FFT begins with the definition of the quantities 


j 


i 


1 
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N , k and r , as follows 

l l 1 


- k - h , -2P” 1 - h 0 *2 P_Z + • • •+ h *2 1 + h »2° 
1 P-1 P _ 2 1 o 


N -I-2P- 1 
1 L 

-2 


and 


r - j 


l 


0 - j_ ,-2 p " 2 +•••+ j »2* + j *2 C 

P -1 2 1 


(4.3.9) 

(4.3.10) 

(4.3.11) 


The exponential term in equation 4.3.8 can be expanded as 
rk (2r ; + j„)ft 4- h p _! 2P- 1 ) 


w N - w 


N 


-W N > W N W 2 


,(4.3.12) 


r h _i 

where the term W 1 p was omitted since it is equal to one. The new 

quantities A are introduced by the definition 

1 

j h , jk 
^ J o P~l J o i 

A (h h ;j ) - Z x(h h ,h )W 2 W N . (4.3.13) 

IP- 2 o o h ^ P-1 l o 

Finally, at the end of the first step of the FFT, equation 4.3.1 
takes the following form 

„ 1 1 1 


r k 
..-* 1 


Z ... r A (h_, h.jj 0 )W 1 .(4.3.14) 


X ( j D _]» • • • » j " XT E t- " '“ n.O 

P 1 ° N h-0h.0h ,.0‘ p2 

o j p-2 

The N values of A can be calculated from the original data directly 
via equation 4.3.13 and require 2N complex operations. The N values of 
X f can be calculated from equation 4.3.14, requiring N»^ complex opera- 
tions. Consequently, the total number of operations is in the first 
step reduced from M 2 to N(^ + 2). 

The former approach can be further generalized. At the step 
where q •* 1,2,.. . ,p-l, the quantities N^, k^, r^ and A^ are introduced 
as follows 


i 


i 


i 


1 


1 



1 


» - Jr!.- 2 p -’ 

q 2 


(4.3.15) 


- k . - h - h „ , -2*~ q ~ l + •••+ h *2 1 + h *2° ,(4.3.16) 

q-1 p-q p-q-1 ! o 


r, ■ - J p _! -ZP-I- 1 +•••+ l q+1 *2 1 + J< 


(4.3.17) 


A q^ h p-q-l * * * * * h o’^q-l » * * * * " 


h -0 


A q-l^ h p-q'* “ ,h o’^q-2 ^ 0 ) w 


Jq -l h p-q Jq-l k q 

W N q-1 . (4.3.18) 


The general expression of the discrete Fourier transform, X r , at 

f-U 

the q step is by induction obtained in the following form 

x <Vi vv ‘ 

r „ k „ 

11 l ^ 

* l z E ••• 2 A q ( Vq-l»“*’ h o :J q-l • (4.3.19) 

h -0 h “0 h -0 4 v 4 4 

c i p-q- 1 

When q is equal to p-1, equation 4.3.19 is reduced to 


r ,k . 

i p- 1 .r-1 

x< Jp-l**** , W " N h f 0 Vl <h o ;j P-2 j o )W Np_1 


. ( 4 . 3 . 20 ) 


The above formulation implies that the discrete Fourier transform, 
X f , of the discrete time history x^, where r,k " 0,1,...,2 P -1 can be 
computed in p steps, from the consecutive calculation of the quanti- 
ties Aq. Each of these steps requires 2N complex multiplications and 
additions, so that the total number of operations is 2pN. Furthermore, 
at each step q, the quantities A can replace the quantities A . in 

q q 

the computer memory, so that only N memory positions are required. 
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: . 3 . 3 A Simple Example 

The following example provides a detailed illustration of the 
FFT approach, as discussed in chapter 4.3.2. A discrete time history, 
consisting of 8 data points, x^, where i * 0,1,..., 7, and extended at 
the time interval [0,T], is given. Following the binary system nota- 
tion, the quantities x^ can be expressed as follows 


x(0,0,0) 

x * x(l,0,0) 


4 

x(0,0,l) 

x - x(l,0,l) 


5 

x (0,1,0) 

x * x(l,l,0) 
6 

x(0, 1 , 1) 

X » x(l,l,l) 


(4.3.21) 

The discrete Fourier transform of this time history consists of 8 
values, X r , where r * 0,1,. ..,7, which, according to equation 4.3.2, 
will have the following expressions in terms of the original data values 


X = -4-(x +x + x +x +x + x +x +x) 

0 8 o 12 3 4 5 6 7 

X = |(x + x W l/8 + x W 2 /® + x W 3 / 0 + x W 4 ^ 8 + x W 5 / 8 + x W 8 / 8 + x W 7 / 8 ) 

1 8 ° 1 2 3 4 5 6 7 

X - T( x + x W 2/ 0 + x W 4 / 0 + x W 6/8 + x + x W 2 / 8 + x W 4 / 0 + x W&/ 8 ) 

2 o ° 1 2 3 45 6 7 

x = + x W 3 ^ 8 + x W^ 8 + X W 1 / 8 + x W 4 / 8 + x W 7 / 8 + x W 2 / 8 + x W 5 / 0 ) 

3 8 O i 2 3 4 5 6 7 

x “ £(x„ + x W 4 / 8 + x + x W 4 / 8 + x + X W 4 / 8 + X + x W 4 / 8 ) 
480! 23 45 67 

X * i(x + x W 5 / 8 + x W 2 / 8 + x W ?/8 + x W 4 / 0 + x W 1 / 8 + x W 6 / 8 + x W 3 / 8 ) 
58 oi 2 3 4 5 6 7 

= |<x 0 + x W 6 / 0 + x W 4 / 0 + x W 2 / 0 + x + x W 6 / 8 + x W 4 / 8 + x W 2 / 8 ) 


X 


6 8 

= + X W 7 / 0 + x W 6 / 8 + x w 5 / 8 + x W 4 / 8 + x W 3 / 8 + x W 2 / 8 + X W 1 / 0 ) 

7 °o 1 2 3 4 5 6 7 

(4.3.22) 

In the above expression, use has been made of the fact that W°, 
where n is an integer, is equal to i. 
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The FFT consists of three steps, l.e. the consecutive calculation 


of , A^, X according to the equation 4.3.18. 
First Step 

j h j k 
1 o 2 2 ? 1 . 

A (h ,h ;j ) - L x(h ,h ,h)W 2 w « 

i i o ° . ~ : • ° 


(4.3.23) 


A (0,0,0) - x(0,0,0)W° + x(l ,0 ,0)W° « x + x 
1 o 4 

A (0,0,1) = x(0,0,0)W° + x(l,0,0)W 4 / 8 - x + x W 4 / 8 
1 ° 4 

A (0,1,0) ■ x(0,0,l)W° + x(l,0,l)W° = x + x 

1 1 5 

A (0,1,1) - x(0,0,l)W l/8 + x(l,0,l)W 5/8 - x W l/8 + x W 5/8 
1 1 5 

A (1,0,0) * x(0,l,0)W° + x(l,l,0)W° ■ x + x 
1 2 6 

A (1,0,1) = x(0,l,0)W 2/8 + x(l,l,0)W 6/8 = x W 2/8 + x W 6/8 

i 2 e 

A (1,1,0) = x(0,l,l)W° + x(l , 1 ,1)W° = x + x 

1 3 7 

A (1,1,1) = x(0, 1 , 1)W 3 / 8 + x (1,1, 1)W 7 / 8 = x W 3 / 8 + x W 7 / 8 (4.3.24) 

1 3 7 

Second Step 

! Vl Va 

a (h ;j ,j ) = I A (h ,h ;j )W 2 W N l (4.3.25) 

2 h^O 1100 

A (0,0,0) = A (0,0,0)W° + A (1,0,0)W° = x +x +x +x 

2 1 1 ° 4 2 6 

A (0,0,1) = A (0,0,1)W° + A (1,0,1)W° = x + x W 4 / 8 + x VJ2/8 + x W 6/8 
2 1 1 4 2 6 

A (0,1,0) = A (0,0,0)W° + A (1,0, 0^/8 = x + x + x W4/8 +• x W 4 /8 


A (0,1,1) = A (0,0,1)W° + A (l.O.DW^/S = x + x W*»/8 + x W6/8 + x W2/8 
2 1 1 4 2 6 

A (1,0,0) = A (0,1,0)^ + A (1,1,0)W° = x + x +x +x 
2 1 1 15 3 7 

A (1,0,1) = A (0,1 ,1)W° + A (1,1,1)W° = x Wl/8 + x W'5/8 + x VJ3/8 + x w?/8 
2 1 1 15 3 7 

A (1,1,0) =■ A (0,1 ,0) VJ2 /8 + A (l,l,0)W6/8 = 

2 1 1 

x W2/8 + x W 2 /8 + x W8/8 + x \]6/B 
15 3 7 

A (1,1,1) - A (0,l,l)W2/8 + A (l,l,l)W6/8 - 
2 1 1 

x W 3 / 8 + x w 7 / 8 + x w 1 / 8 + X w 5 / 8 (4.3.26) 

1 5 3 7 


(4.3.26) 
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Third Step 


X(j ,j ,j ) - T EA(h:J ,j )W K 2 
2 l o 8 v 2 0 1 o 


(4.3.27) 


X - X(0,0,0) - -4A (0,0,0)W° + A (1,0,0)W°] - 

0 o 2 2 

T 

“ «( x + x + x + x + X v + X + x + x ) 

8 ° 4 2 6 1 5 3 l' 

X, - X(0,0,1) - -[A (0,0,1)W° + A (1,0, 1)W ] - 

1 8 2 2 

- J(x n + x W 4/8 + x W 2 / 8 + x W 6 / 8 + x W 1/8 + x W 5/8 + x W 3 / 8 + x W 7 / 8 ) 

8 o 4 2 6 1 5 3 7 

X, = X(0,1,0) * |[A (0 ,1 ,0)W° + A (1,1,0)W°] - 

i 0 2 2 

= £(x + x + x W 4 / 8 + x W 4 / 8 + x W 2 / 8 + X W 2 / 8 + x W 6/8 + x W 6 / 8 ) 

8 o 4 2 6 1 5 3 7 

X, = X(0,1,1) = 1-lA (0,1,1)W° + A (1,1,1)W°] - 

3 a 2 2 

= 7 (x + x W 4 / 8 + x W 8 / 8 + x W 2 / 8 + x W 3 / 8 + x W 7 / 8 + X W 1 / 8 + x W 5 / 8 ) 

8 0 4 2 6 1 5 3 7 

X, = X(1,0,0) = TfA (0,0,0)W° + A (1,0,0)W 4/8 1 = 

•+ « 2 2 

= ^(x + x + x + x +x W 4 / 8 + x W 4 / 8 + x W 4 / 8 4 x W 4 / 8 ) 

8 o 4 2 6 1 5 3 7 

X. = X(1,0,1) = J[A (0,0,1)W° + A (1,0,1)W 4/8 ] = 

D 0 2 2 

= i(x + x W 4/8 + x W 2/8 + x W 6/8 + x W 5/S + x W l/8 + x W 7/8 + x W 3 ^ 8 ) 
o O 4 2 6 1 5 3 7 

X 6 = X(1,1,0) = |{A 2 (0,l,0)W O + A 2 (1,1,0)W 4/s ] * 

= f(x + x + x W 4 ^ 8 + x W 4 ^ 8 + x W 6/8 + x W 6 ' 8 + x W 2 ' 8 + x W 2 / 0 ) 
8 042 6 1 5 3 7 

X 7 = X(l,l,l) = J[A (0,1,1)W° + A (1,1,1)W 4/8 ] - 

/ o 2 2 

= £(x + X W 4 ' 8 + x W G/,S + x W 2 ' 8 + x W 7 / 8 + x W 3/8 + x W 5 / 8 + X U 1 ^ 8 ) 

° O 4 2 6 1 5 3 7 

(4.3.28) 

The final expressions In the equation 4.3.28 are identical to the 
expressions in equation 4.3.22. Nevertheless, the required number of 
operations with the FFT method is 24, while, according to equation 4.3.22, 
the number of operations would be 64. The order of operations and the 
elements used at each step of the previous example are illustrated in 
figure 15. 
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4.3.4 The FFT Computation 


Although a FFT subroutine was available from the V.P.I. & S.U. 
computer library, it proved most efficient to develop a new program 
based on the formulation of chapter 4.3.2. The flow chart of this 
program is presented in figure 16. 

A considerable saving in the computer time was obtained by the 
recursive calculation of the sines and cosines, in which the expo- 
nential terms are expanded, instead of recalling the built-in sub- 
routines. The general exponential term is expanded as 
rk rk 

N n 2Trrk 2irrk 

W =* e = cos — i*sin — — , rk = 0,1,. .. ,N-1 , (4.3.29) 


N 


N 


while the sines and cosines can also be expanded as 


2nrk 

cos — — * cos 


27r(rk-l) 2ir 2ir(rk-l) 2ir 

N cos “N ’ sin N sin "i 


and 


2irrk 2n(rk-l) 2ir 2ir(rk-l) 2ir 

sin — jj” - cos N sin N + sin N cos — 


(4.3.30) 


According to the expressions 4.3.30, all sines and cosines 
required can be calculated with simple operations from the initial 

2tt 2 

values of cos and sin "m. Furthermore, since all sines and cosines 
N " 

are simply related to the sines and cosines in the first quadrant, 
the expressions 4.3.30 need to be used only for values of rk equal to 

- i 


m 


-1 — 


„wL 
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4.4 Spectral Calculations In Each Block of Data 


The FFT technique, as discussed in chapter 4.3, was developed for 
complex data sequences. When applied to real time histories, such as 
the records of atmospheric measurements, it permits the simultaneous 
computation of the Fourier transforms of two of these real time his- 
tories. 

Given the two real discrete time histories, x^ and y^» where 
k * 0,1,..., N-l, a new complex sequence, z^, can be generated as 
follows 


Xj^ + iy fc , k - 0,1,..., N-l 


(4.4.1) 


If Z , r ■ 0,1,..., N-l, is the discrete Fourier transform of z 

V 

and Z* its complex conjugate, it has been proved (see for instance 
reference 1), that the discrete Fourier transforms of the time his- 
tories x r and y r can be respectively calculated as 

1 n 

X r * 2 [Z r + Z N-r ] » r-O.l,..,,^ (4.4.2) 


T r • 2? Z r - ’ r ' °» l f 


(4.4.3) 


It has to be mentioned that the finite Fourier transform of a 
real time history has independent values only in the first half of 
the frequency domain, namely for frequencies up to the so called 
Nyquist frequency — . The values of the Fourier transform in the 
second half of the frequency domain are simply the complex conjugates 
of respective values in the first half, i.e. 
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X r - X*_ r , r - |,| + 1 N . (4.4.4) 

Taking advantage of the relations 4.4.2 and 4.4.3, the discrete 
Fourier transforms of two of the quantities T, u, v and w at a time 
are simultaneously calculated for each block of 8192 data points. The 
block-spectral density function estimates will be calculated with use 
of the expressions 4.1.1 and 4.1.4, appropriately modified for discrete 
time histories. 

Let X n and Y^, r ■ 0,1,...™ , be the discrete Fourier transforms 
r r 1 

of the n block of the blocked time histories x^ and y^, k * 0,1,..., N-l. 
The block-power spectral density function estimate, G^(f r ), is then 
calculated as 

G x (f r> " fl X ?l 2 * (4 - 4<5) 

and the block-cross-spectral density function estimate, as 

G xy Cf r ) " f X r* Y r » (4.4.6) 

N 

where r - 0,1, ... and n - 1,2,... ,M. 

The estimates of the coincidence and quadrature spectral density 
functions for each block can be independently calculated as 

C |J tfrJ »f[Re(x“) Re(Yj) - Im(x“) Im(Yj)] (4.4.7) 

and 

-Q£y( f r ) = f(Re(x") I«(Y“) + Im(xJ) Re(Y^) ] . (4.4.8) 

All block-spectral estimates need further to be multiplied by 
the scale factor 1/0.875, due to the cosine taper data window. 
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4.5 Smoothing of the Spectral Estimates 

n n 

The block-spectral estimates, G x and G^, as calculated in chapter 

4.4, follow chi-square distributions of 2 degrees of freedom with 

mean values the respective spectral density functions, G or G , and 

a xy 

standard deviations equal to the respective mean values (see reference 
1) . This implies that the average systematic error of the estimate 
will be equal to 100Z, which is quite unacceptable. Consistent esti- 
mates of the spectral density function can be obtained by smoothing 
the block estimates. A combined smoothing technique, of first smooth- 
ing over the ensemble of estimates (or segment averaging) and the 
smoothing over appropriate frequency intervals will be used to pro- 
duce appreciably accurate spectral estimates. The combined smooth- 
ing is performed for the power spectral estimates and for the co- and 
quad-spectral estimates, which are of main interest. 

Let G represent any of the quantities G x , Cjjy and Q^. The en- 
semble smoothing consists in averaging all block spectral estimates, 

G n , so that the respective sample spectral estimates, G, are obtained, 
.e. 


1 M n 

G(f ) - m E G < f r > 

r n»l 


r - 0,1 j " 1 


(4.5.1) 


The frequency smoothing consists in averaging all sample spectral 

estimates over a frequency interval [f^^.f^^] and assign the average 

values to the center frequency f . . Therefore, each value of the 

„ r ^2 

smoothed spectral estimates, G, stands for l values of the raw spec- 
tral estimates, G, and is generated according to the expression 
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1 

£ 


m*l 


(4.5.2) 


The nonuniform distribution of energy in the spectra requires 
that the frequency Intervals for frequency smoothing are longer at 
the high frequency regions, where less power is present. The pro- 
posed frequency smoothing provides 46 values of the smoothed spec- 
tral estimates out of 4096 values of the raw spectral estimates. The 
calculation of these values and their distribution in the spectrum 
are obtained from the following expressions 


q « 1.....4 , £ - 1 

G(f r ) = G(f r ) 
q * 5 , ... ,15 , £ * 4 


r - q 


, (4.5.3) 


1 ^ — 

G(f_j_ 0 ,) - 7 Z GU^) 


r+2.5 


m*l 


q * 16 .... ,24 , £ - 16 


16 


g <W 5 > ‘ * E , G <W 

m»=l 


q ■ 25 33 , £ * 64 

1 64 


, r ■ 4(q-4) , (4.5.4) 


, r = 16(q-13) , (4.5.5) 


G(f r + 32.5 ) " M * r ’ 64( «‘ 22 ' * (4 - 5 ‘ 6 > 


q * 34, . . . ,46 , £ * 256 


256 


G<£ rtl28.5 ) - 256 * r ‘ 256( ’- 31) • < 4 - 5 ' 7) 

m-x 

As discussed before, the estimates of the smoothed power, G x , 
coincidence, C x y, and quadrature, Q x y, spectral density functions are 
directly calculated with the use of the expressions 4.5. 1 , . . . ,4.5. 7. 

The smoothed estimates of the remaining spectral quantities, as defined 


in chapter 4.1, have still to be computed. The magnitude, |G X y(f) |, 
and the phase angle, 0 X y(£), of the smoothed cross-spectral density 


estimate are respectively calculated as 


<V f >! " <0 


(4.5.8) 


e xy (f) - arctan 


C xy (f) 


(4.5.9) 


while the smoothed estimate, Y^y(f), of the coherency function is 


calculated as 


^y< f) 


|G xv (t )| 2 


G x (f) Gy(f) 


(4.5.10) 


4.6 Plotting of the Spectra 


In this report, the product of the spectral density function at 
a certain frequency and the corresponding frequency itself is plotted 
versus the natural logarithm of the frequency. This type of presenta- 
tion of the spectra provides a more detailed picture of the spectral 
functions in the higher frequency region where much less power is 
present with comparison to the low frequency region. In addition, this 
way of spectral plotting provides for a method of fast checking of the 
spectral estimates. 

The variances and covariances of a pair of mean-free time histories 
are related to the power and coincidence spectral density functions 
respectively (see reference 1) by the relations 


x* - / G (f ) df 
0 x 


(4.6.1) 
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and 


xy - / C (f)df . (4.6.2) 

0 J 


By simply transforming the integrand quantities, one may rewrite 
the above relations as 


** - / f *G (f) d(lnf) (4.6.3) 

0 


and 


xy - / f*C xV (f) d(lnf) . (4.6.4) 

o y 


Therefore, the accuracy of the spectral estimates can be checked by 
comparing the respective values of the variances and covariances, as 
calculated with procedures described in chapter III, to the areas under 
the curves obtained by plotting the product of the frequency and the 
power or coincidence spectral density function values at this frequency 
versus the natural logarithm of this frequency. 
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CHAPTER V 

DISCUSSION OF THE RESULTS 

The objective of this thesis is to provide an efficient method 
for the measurement and the analysis of atmospheric turbulence data 
and specifically to explain the procedures used for the statistical 
and spectral analysis of the data. Therefore, the discussion is based 
on one complete set of results obtained from the analysis of the data 
recorded at one specific run with one of the anemometer svstems. From 
the several records of data which were available at the time when this 
thesis was written the data of run number 1 were selected. These data 
were obtained on August 20, 1973, with the TSI #1193 anemometer at the 
10C- foot level of the 250-foot meteorological tower located at NASA 
Wallops Island Flight Center. Calibration specifications for the TSI 
#1193 anemometer are provided in figures 7, 8, 9, 10 and 11 and in 
table I. For the mean velocity and mean probe-yaw angle of run 1, 
shown in table IV, it can be seen from the figures 11 and 10 respec- 
tively that the expected error in the calculation of the magnitude 
of the velocity is of the order of 1 percent while the expected error 
in th<" calculation of the velocity components does not exceed 4 percent. 

The data-record consists of 475,136 data points, recorded at a 
time increment of 0.005 seconds over a total time interval of approxi- 
mately 39 minutes and 35 seconds. As discussed in section 3.3, it is 
divided into 58 blocks consisting of 8192 data points each, extending 
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over a time interval of approximately 41 seconds. The block-mean values 
for the velocity components in the sensor oriented coordinate system, 
their standard deviations, the velocity magnitude and the probe-yaw 
angle are shown in table II. By inspecting the variation of these quan- 
tities over the total number of blocks, one may observe a random rather 
than systematic varying distribution. Only for the last four blocks 
the mean wind seems to change slightly in magnitude as well as direction 
but without any considerable change in its frequency composition. The 
total number of reverse arrangements for the block-means of the above 
quantities is also shewn in table II. One can see that all the calcu- 
lated values for the number of reverse arrangements lie within the limits 
of assumed stationarity as discussed in section 3.4. Therefore, it 
was decided to consider the data of run 1 as being stationary and, hence, 
to proceed with the statistical and spectral analysis. 

The block-means of the temperature and the velocity components in 
the mean-wind oriented coordinate system are shown in table III. It can 
be seen from this table that the vertical and lateral components take 
values much smaller than the values of the longitudinal component, which 
confirms that the mean wind is nearly horizontal, and nearly parallel to 
the axis of the probe. The sample-mean values of the temperature, the 
velocity components in the mean wind direction, the magnitude and phase 
angle of the velocity vector and all variances and covariances, obtained 
by averaging the respective block-mean values, are shown in table IV. 

As shown in this table, the calculated value of the sample probe-yaw 
angle is very small (approximately 1.34 degrees), which is a good indi- 
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cation for the accuracy of the measurements. This mean wind is coming 
from a south-eastern direction, namely di-ectly f rr>m the Atlantic Ocean 
without any considerable upstream obstacle distortion. 

It has to be emphasized at this point that the values of the vari- 
ances and covariances depend on the length of the block-time interval 
considered, due to the change in the block-mean values which are sub- 
tracted from the instantaneous values. As discussed in section 4.6, 
the variances and covariances represent the area under a spectral curve. 
This area obviously depends on the frequency interval considered. For 
instance, it has been generally accepted that the complete power spectrum 
of the longitudinal wind component has a form similar to the one shown 
in figure 17, as calculated by Van der Hoven from measurements at Brook- 
haven. This spectrum extends over a wide range of frequencies and pre- 
sents two major peaks. The peak at the low-frequency region (at a period 
of approximately 100 hours) corresponds to periodic passage of cold or 
warm fronts, while the peak at the high-frequency region (at a period of 
the order of 1 minute) corresponds to the turbulent fluctuations of the 
wind. It Is the latter part of the velocity spectrum, often called the 
"gust velocity" spectrum, which is of incerest in the study of turbulence. 
The region of the spectrum between the two peaks, called the "mesometeor- 
ological gap" contains very little power, due to the fact that no normal 
physical process can produce velocity fluctuations with frequencies 
lying in this range. 

Therefore, the values of variances and covariances are consistent 
and romparable to each other, only when they are calculated from data- 
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blocks extending over time intervals large enough to determine the 
spectrum down to frequencies near the mesometeorological gap, usually 
of the order of one cycle per hour. This, of course, would require 
data-saraples extending over intervals of the order of several hours. 

For a given fixed data-sample, the block-time length can be increased 
by dividing the sample into a smaller number of blocks. This method 
has the disadvantage that it provides less accurate values for the 
sample-variances and covariances, since they are calculated as the 
average of a smaller number of block-values. 

The block-time length effect can be well illustrated with the 
use of the turbulence intensities, defined as the ratios of the stan- 
dard deviations of the velocity components in the mean wind oriented 
coordinate system and the velocity magnitude. The turbulence infen- 
sities calculated from the data of run 1 for the block time lengths 
of approximately 41, 328 and 1188 seconds are shown in table V. It 
can be seen from this table that the values of the turbulence inten- 
sities increase considerably as the block-time length is increased, 
especially in the case of the lateral velocity component. 

The block-time length is also of great importance for the power 
and cross-spectra, since it determines the lower frequency for which a 
spectral density fu> .cion can be estimated from a given sample of data. 
Consider a sample of data consisting of M blocks of N data points each, 
sampled at a time increment At. Then the block-time length is equal to 
N*At, while the total sample-time length will be equal to M«N*At. The 
smooth spectral estimates, as discussed in chapter IV, can be obtained 
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for the frequency range between the lowest frequency (or fundamental 
frequency), 1/N.At and the Nyquist frequency (or folding frequency), 

1/2* At. Since the number of data per block, N, Is held fixed equal 
to 8192 due to the FFT needs, the only way of estimating the spectra 
in lower frequency regions is to increase the time increment. At, i.e. 
to consider selectively only one part of the data points, separated from 
each other with a longer time increment. At. As a result of the fixed 
sample-time length, the number of blocks, M, has to be adequately de- 
creased, with the consequence of decreasing the effectiveness of the 
ensemble smoothing. 

The smooth spectral estimates in the frequency range from approxi- 
mately 0.0244 hertz up to 100 hertz, are calculated by considering all 
data points, namely taking At = 0.005 seconds and M ■ 58. Spectral 
estimates for lower frequencies down to 0.0030 hertz are obtained by 
considering one data-point out of eight, i.e. by taking At * 0.040 
seconds and, as a consequence, M = 7. An attempt of estimating spec- 
tral values for frequencies as low as 0.00084 hertz by taking At = 0.145 
seconds and M = 2 (i.e. considering only one every 29 data points) was 
abandoned, as producing unacceptably high distorted spectral estimates 
as a result of the poor ensemble smoothing. 

The power spectra of the longitudinal, the lateral and the verti- 
cal velocity components and the temperature are shown in figures 18, 

19, 20 and 21 respectively. It can be seen in these figures that the 
power spectral estimates show a considerably small amount of scatter, 
expecially for frequencies higher than approximately 0.1 hertz. At the 
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lower frequency region, there exists higher scatter of the spectral 
values, due primarily to the absence or Inefficiency of the frequency 
smoothing, as discussed In section 4.5. It can also be seen that spec- 
tral values estimated from data sampled at time Increments equal to 
0.005 and 0.040 seconds are In general consistent, while Increasing the 
time increment to 0.145 seconds provides highly distorted values. The 
form of the curves that can be fitted to the calculated spectral values 
for the velocity components is in general monotonously decreasing, which 
confirms that the estimated spectra lie in the gust velocity region. The 
spectral curves approach zero at some frequency higher than 100 hertz. 
Some high power spectral values, not shown in the above figures, were 
calculated for frequencies close to 100 hertz but were discarded as 
resulting from resolution errors. It has to be recalled at this point 
that all frequencies higher than 100 hertz were removed from the original 
data with the use of sharp low-pass electronic filters. 

The estimated values of the magnitude of the cross-spectral density 
function, the co-spectral density function and the coherence function 
between the longitudinal and the vertical velocity components are shown 
in figure^ 22, 23 and 24 respectively. These estimates are in general 
less smooth than the power spectral estimates, due to the higher eri r 
involved in the calculations between two time histories. Of specific 
interest is figure 24, since it indicates that the longitudinal and the 
vertical velocity component fluctuations are completely incoherent for 
all frequencies below a cut-off frequency of approximately 5 hertz. 
Nevertheless, these fluctuations are considerably coherent for higher 
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frequencies, with an average value of the coherence function approxi- 
• mately equal to 0.6. 

The checking of the accuracy of the spectral estimates, as sug- 
gested in section 4.6, was performed for all the above cases with 
| satisfactory results. 




CHAPTER VI 

CONCLUSION AND RECOMMENDATIONS 

Concluding Lite discussion that has been made in the previous chap- 
ters of this thesis, one could summarize what was achieved as follows. 

Highly accurate measurements are necessary for the detailed des- 
cription of the wind flow in the atmospheric boundary layer. The data 
obtained from these measurements have also to be in an appropriate form 
in order to provide accurate estimates of the required statistical and 
spectral quantities. To meet these requirements, a complete instrumen- 
tation system was designed and a new data-analysis method was developed. 

The TSI 1080 three-dimensional split-rilm anemometer system was 
selected for the measurement of the instantaneous values of the wind- 
velocity vector and the temperature. It provides seven output voltages 
from which the velocity components in the sensor-oriented coordinate 
system and the temperature can be calculated. The general calibration 
procedure, which was suggested by TSI, was proved to be complicated 
and inaccurate and, therefore, a new co e calibration procedure had 
to be developed. This calibration was based on the recognized fact that 
the anemometer gives the most accurate results when the probe operates 
with its axis parallel to the direction of the mean-wind vector and 
the sensors directed into the wind. 

All measurements that were necessary for the calibration were per- 
formed at the Quality Verification and Calibration Facility at NASA 
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Wallops Flight Center. The entire procedure was streamlined Into a 
standard form, allowing a quick and efficient calibration of the large 
number of anemometers used In this research program. Based on this 
calibration, a computer program was developed, which converts the out- 
put voltages of the anemometers into the actual velocity components 
and the temperature. The accuracy of the obtained quantities was 
checked and found satisfactory. 

Six probes were mounted at the elevations of 30, 50, 100, 150, 

200 and 250 feet on the 250-foot meteorological tower located at Wallops 
Island on the Atlantic coasc. All probes could be rotated so that their 
axis could be placed closely parallel to the mean wind direction. A 
new complete data-acquisition system was designed and fixed in an instru- 
mentation trailer, which enabled easy transportation and protection of 
the system. At each run of the system, the output voltages from all 
six anemometers were simultaneously recorded on analog magnetic tapes. 
Later, each analog tape was played back and the analog signals from each 
anemometer were digitized, sampled at a rate of one data-point every 
0.005 seconds and recorded on a digital magnetic tape. Due to its 
complexity, the data-acquisition system requires elaborate calibration 
checks to insure its proper operation* The time duration of each data- 
sample is approximately 40 minutes, as imposed by the maximum capacity 
of the analog tapes. It has been attempted to increase the sample time 
by using a lower speed of the analog tape-recorder, but this resulted 
in unacceptably high noise-levels in the recorded signals. 

Each of the above digital magnetic tapes contains a record of the 
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seven output voltages of one TSI anemometer at one run of the system. 
With the use of the program "CONVERSION", as shown In the appendix, 
a new magnetic tape was generated, containing a record of the three 
velocity components in the sensor oriented coordinate system and the 
temperature. Each time history consists of approximately 480,000 data 
points, separated from each other with a time interval of 0.005 seconds. 

Consistent estimates of the statistical and spectral quantities 
could be obtained only if certain averaging procedures would be applied 
to the data. Therefore, each time history was divided into a number 
of data-b locks, so that all desired quantities would be calculated as 
averages of the respective block-values. The number of data-points in 
each block was fixed and equal to 8192, as imposed by che fast Fourier 
transform needs and the capacity of the available computer. The entire 
analysis was based on the assumption that the measured time histories 
were stationary, which implies that all statistical quantities do not 
vary significantly over the total sample time interval. This assump- 
tion was checked with the use of the program "STATIONARITY" . This 
program provides the block -mean values and the block-standard deviations 
of the velocity components and the temperature as well as the results 
of a statistical test which detects monotonic trends of the above quan- 
tities. Each data-record was further processed, only if the results of 
the inspection of the statistical block-values and the trend test would 
ensure stationarity within reasonable limits of confidence. 

The velocity components in the sensor-oriented coordinate system 
do not possess any physical meaning. Therefore, a coordinate transfor- 
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mation Into a coordinate system connected with the mean wind direction 
was necessary. This transformation provided the longitudinal, lateral 
and vertical velocity components, which could give a good description 
of the properties of the wind velocities. The new velocity components 
were calculated and together with the temperature recorded on a new 
digital magnetic tape with the use of the computer program "MEAN-WIND". 
This program also provides the variances and covariances of the above 
four time histories as well as the turbulence intensities. As already 
mentioned, each time history has beer divided into a number cf blocks 
and all sample-quantities were calculated as the averages of the block- 
values . 

One of the main objectives of this thesis was to provide a method 
for the calculation of reasonably accurate measurements of the power 
and cross-spectral density functions of all three velocity components 
and the temperature, expecially in the high frequency range between 
0.1 and 100 hertz. Such a high frequency range is necessary for the 
understanding of the turbulent stresses, turbulent transport phenomena 
and turbulent dissipation. Spectral estimates for lower frequencies 
down to 0.003 hertz were also calculated, but with lower accuracy due 
to the limited time length of each data-sample. 

It was recognized that the finite record-length of the data used, 
would produce a certain amount of distortion in the spectral estimates. 
Therefore, the cosine taper data window was used to reduce this effect. 
The power, cross-, co- and quad-spectral density function estimates and 
the coherence function estimate were calculated with the direct Fourier 
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transform method. A new fast Fourier transform technique was developed 
and used to enable fast and effective computations. This FFT technique 
possesses the addidional advantage of not requiring any bit-reversal 
operation. A recursive computation of the necessary sines and cosines 
was also used to further reduce the computation time. 

First, a set of raw spectral estimates were obtained for each 
block of 8192 data-points. Each block-spectral estimate consists of 
4096 discrete values. Then all block-spectral values corresponding 
to each frequency were averaged to produce the sample-spectral estimates, 
consisting also of 4096 discrete values. Finally, the sample-spectral 
values were averaged over appropriate frequency intervals to produce 
the smooth spectral function estimates, consisting of 46 spectral values 
each. The frequency intervals were arranged in such a way, that the 
consequent smooth spectral values would be separated by an. approximately 
constant distance, when plotted versus a logarithmic frequency scale. 

It was also decided to represent the spectral estimates by plotting the 
smooth spectral density functions multiplied by the respective frequency 
versus the natural logarithm of the frequency. This way, a clear repre- 
sentation of the spectra in the high frequency region was obtained, 
although very small power is present in this region. A fast way of 
checking the accuracy of the spectral curves was possible by comparing 
the areas under these curves to the respective values of variances and 
covariances. The smooth spectral estimates are in general much more 
reliable than the raw spectral estimates, since they are almost free of 
scatter which may occur in a specific sample but is not due to an 
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existing physical process. The smooth estimates obtained according to 
the methods discussed, are considered to be reliable for use in most 
applications. All computations that were necessary in order to obtain 
the smooth spectral estimates of the three velocity components and the 
temperature, were performed with the use of the program "SPECTRA", also 
shown in the appendix. This program can be used in any sufficiently 
large modern digital computer. 

The discussion in this thesis was confined on one single record 
of data, obtained with one anemometer system at the 100-foot level of 
the meteorological tower. The main concern was to develop and check 
the proposed method for the measurement and the statistical analysis 
of atmospheric turbulence. As soon as several complete sets of measure- 
ments from all six anemometers at the different elevations of the 
meteorological tower are available, this method can be used to provide 
sufficient material for a documented discussion on several important 
aspects of the turbulence in the atmospheric boundary layer. 

At the time this thesis was written, a number of modifications and 
extensions of the above method were conceived but were not executed due 
to the lack of time and the limitations of the computing facilities at 
NASA Wallops Flight Center. Nevertheless, these modifications and 
extensions are presented here so that they can be used in the continu- 
ation of this research project. 

1. Increase of the number of data per block . A final modification 
of the computer program used for spectral calculations enabled the reduc- 
tion of the amount of computer memory required, so that a number of 
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16,384 data points could be simultaneously handled with the IBM 370 
computer. Therefore, the number of data per block can be increased from 
8,192 * 2 13 to 16,384 * 2 1 **. This modification will increase the fre- 
quency interval for which spectral estimates can be obtained to the 
frequency range between approximately 0.0122 and 100 hertz and is also 
expected to improve the accuracy of all spectral estimates. 

Calculation of the auto- and cross-correlation functions . As 
discussed in section 4.1, the auto-correlation function and the cross- 
correlation function form Fourier transform pairs with the pcwer and the 
cross-spectral density functions respectively. Therefore, the correla- 
tion functions, R^t) and R^ix), can be calculated as the inverse 
Fourier transforms of the spectral density function estimates, G (f) 

A 

A 

and G xy (f), respectively. 

When auto- and cross-correlation functions are calculated from 
spectral estimates, obtained from finite discrete time histories, cer- 
tain considerations have to be made in order to avoid erroneous results. 
Details on these considerations and computational recipes can be found 
in the ""eferences 1 and 16. 

3. Calculation of the integral scales of turbulence . The integral 
scale lengths of turbulence characterize the maximum distance of two 
points whose velocity fluctuations are correlated. Once the auto-cor- 
relation functions are calculated, as discussed in the previous section, 
simple integrations can provide the Eulerian integral time scales of 
turbulence in the mean wind direction. Then, by assuming the validity 
of Taylor's hypothesis, the integral scale lengths of turbulence in 
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the mean wind direction can also be obtained. These scales are an 
estimate of the physical size of the turbulent eddies in the direction 
of the mean wind at the different tower levels. 

4. Calculation of the cross-spectra and cross-correlations between 
the velocity components and the temperature at different altitudes . 

These quantities can be easily calculated by combining the data obtained 
from the anemometers located at different elevations on the meteorological 
tower. The study of the variation of these quantities with height can 
provide significant information about the vertical momentum and heat 
flux as well as the vertical scales of the turbulence. 
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Figure 4: Multiplexing and Analog Recording System. 
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Figure 5: Demultiplexing, Digitizing and Digital Recording System. 
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TABLE II. BLOCK-MEAN VALUES OF THE VELOCITY COMPONENTS III THE SENSOR ORIENTED COORDINATE SYSTEM, THEIR 
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TABLE III. BLOCK-MEAN VALUES OF THE VELOCITY COMPONENTS 
IN THE MEAN WIND ORIENTED COORDINATE SYSTEM 
AND THE TEMPERATURE (RUN 1, lOO-FOOT LEVEL) 


Block 

Number 

Temperature 

Velocity Components 

n 

tVf) 

u^fps) v^fps) w^fps) 


1 

82.14 

14.34 

0.75 

0.40 

2 

85.22 

13.61 

0.26 

0.78 

3 

85.67 

13.08 

1.09 

0.47 

4 

83.44 

13.50 

0.45 

0.49 

5 

80.73 

13.32 

-1.72 

1.70 

6 

81.67 

13.69 

0.31 

0.38 

7 

82.20 

14.82 

-0.33 

0.98 

8 

85.88 

13.66 

0.24 

0.88 

9 

80.65 

13.65 

0.64 

0.25 

10 

83.21 

13.64 

-0.58 

0.17 

11 

81.99 

14.74 

-0.98 

0.86 

12 

82.76 

14.68 

-1.62 

1.14 

13 

81.99 

13.82 

-2.67 

0.18 

14 

81.43 

14.41 

-2.30 

0.51 

15 

81.24 

14.04 

-1.45 

0.23 

16 

82.08 

13.63 

-1.07 

0.49 

17 

81.20 

13.99 

-1.95 

0.18 

18 

81.36 

14.48 

-1.49 

0.82 

19 

81.48 

14.45 

-1.24 

0.80 

20 

82.39 

13.22 

-1.55 

0.25 

21 

81.23 

13.49 

-0.85 

0.06 

22 

83.07 

13.26 

0.29 

0.19 

23 

89.02 

13.34 

0.50 

0.72 

24 

84.26 

13.66 

0.35 

0.41 

25 

82.53 

13.52 

-1.64 

0.98 

26 

82.27 

13.44 

-1.78 

0.60 

27 

79.82 

13.14 

1.27 

0.38 

28 

80.50 

12.28 

2.03 

0.08 

29 

86.61 

12.93 

1.57 

0.77 

30 

82.65 

13.92 

1.51 

0.01 

31 

85.08 

14.37 

1.64 

0.57 

32 

86.93 

14.17 

1.56 

1.06 

33 

84.41 

14. Cl 

0.94 

0.52 

34 

84.43 

14.17 

-0.40 

0.89 

35 

84.31 

14.29 

0.05 

1.01 

36 

84.46 

13.74 

-0.14 

1.01 

3/ 

85.67 

13.36 

0.34 

0.11 

38 

81.97 

13.65 

- 0.88 

0.81 



!? 



TABLE III. (Continued) 


i 


Block 

Number 

Temperature 

Velocity Components 

n 

t“( # p) 

u^fps) v^Cfps) w^Cfps) 


39 

81.15 

12.74 

-0.87 

-1.33 

40 

83.32 

12.61 

-1.73 

-0.97 

41 

81.78 

13.75 

-1.42 

-0.72 

42 

81.06 

12.84 

-1.63 

-0.81 

43 

82.08 

12.54 

-1.12 

-0.84 

44 

81.58 

12.70 

-0.97 

-0.84 

45 

81.92 

13.34 

-1.21 

-0.88 

46 

84.78 

12.37 

-3.10 

-1.01 

47 

81.80 

12.20 

-1.21 

-0.60 

48 

81.22 

13.51 

-1.95 

-1.02 

49 

82.57 

14.30 

C .64 

-0.04 

50 

80.22 

14.56 

0.75 

0.46 

51 

81.03 

13.98 

-0.32 

0.25 

52 

84.12 

15.79 

-0.14 

-1.05 

53 

84.07 

15.36 

-0.16 

-0.69 

54 

81.28 

15.88 

2.03 

-0.30 

55 

77.92 

15.80 

4.66 

0.23 

56 

77.31 

16.36 

6.16 

-0.25 

57 

77.20 

16.49 

5.83 

-0.49 

58 

77.33 

16.61 

4.58 

0.06 
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TABLE IV. SAMPLE MEANS, VARIANCES AND COVARIANCES OF THE TEMPERATURE 
AND THE VELOCITY COMPONENTS IN THE MEAN WIND ORIENTED 
COORDINATE SYSTEM (58 BLOCKS, RUN 1, lOO-FOOT LEVEL) 


Sample Means Sample Variances Sample Covariances 


u ■ 

13.90 

fps 

u 7 * 

0.475 

(fps) 2 

tiv 

- 

0.063 

(fps) 2 

V = 

0.00 

fps 


0.753 

(fps) 2 

uw 

ss 

-0.102 

(fps) 2 

w = 

-0.53 

fps 

^ = 

0.391 

(fps) 2 

u6 

- 

.0.178 

fps • °R 

U = 

13.91 

fps 

Q 7 - 

6.46 

(°R) 2 

vw 

3 

0.060 

(fps) 2 

S ■ 

1.34 

degrees 




v6 

- 

-0.111 

fps*°R 

f = 

82.44 

°F 




w9 

a 

-0.390 

fps * °R 







1 1 

1 i 
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TABLE V. 

TURBULENCE INTENSITIES 
(RUN 1, lOO-FOOT LEVEL) 

FOR DIFFERENT BLOCK 

TIME LENGTHS 

« 

T(sec) 

p/u 

p/u 

fF/U 


*41 

0.0495 

0.0624 

0.0450 


~328 

0.0648 

0.1101 

0.0527 


-1188 

0.0871 

0.1455 

0.0556 



APPENDIX 

LISTING OF FORTRAN PROGRAMS 

I. CONVERSION 

II. STATIONARITY 

III. MEAN-W7.ND 


IV. SPECTRA 
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